#include "kernel/solution.hpp"

Component:: Component( const Library &library, const string &name ) :
C(0),
species( library[name] )
{
	if( !library.is_closed() )
		throw exception("Component('%s'): Library is NOT closed", name.c_str() );
	
}


Component:: ~Component() throw()
{
}

Component:: Component( const Component &other ) throw() :
C( other.C ),
species( other.species )
{
	
}

const string & Component:: key() const throw() { return species.name; }

Solution:: ~Solution() throw()
{
}

Solution:: Solution( const Library &lib ) :
components_( lib.size(), as_capacity ),
size( lib.size() ),
library( lib )
{
	if( !library.is_closed() )
		throw exception("Solution(Library is not closed!)");
	
	for( size_t i=1; i <= size; ++i )
	{
		const Species &s = library(i);
		const Component comp( library, s.name );
		if( !components_.insert(comp) )
			throw exception("Unexpected multiple component '%s'",  s.name.c_str() );
	}
	
	
}

Real & Solution:: operator[]( const string &id ) 
{
	Component *s = components_.search( id );
	if( !s ) throw exception("no Solution['%s']", id.c_str() );
	return s->C;
}


const Real & Solution:: operator[]( const string &id ) const
{
	const Component *s = components_.search( id );
	if( !s ) throw exception("no const Solution['%s']", id.c_str() );
	return s->C;
}

Real & Solution:: operator[]( size_t index ) throw()
{
	assert(index>0);
	assert(index<=size);
	return components_(index).C;
}

const Real & Solution:: operator[]( size_t index ) const throw()
{
	assert(index>0);
	assert(index<=size);
	return components_(index).C;
}


const Species & Solution:: operator()(const string &id ) const
{
	const Component *s = components_.search( id );
	if( !s ) throw exception("no Solution('%s')", id.c_str() );
	return s->species;
}

const Species & Solution:: operator()( size_t index ) const throw()
{
	assert(index>0);
	assert(index<=size);
	return components_(index).species;
}


void Solution:: load_from( const ChemSys &chemsys )
{
	if( chemsys.M != size || &library != &chemsys.library )
	{
		throw exception("Solution.load_from(mismatching chemsys)");
	}
	ComponentSet::iterator c = components_.begin();
	for( size_t i=1; i <= size; ++i, ++c)
	{
		(*c).C  = chemsys.X[i];
	}
}

void Solution:: save_into( ChemSys &chemsys ) const
{
	if( chemsys.M != size || &library != &chemsys.library )
	{
		throw exception("Solution.load_from(mismatching chemsys)");
	}
	ComponentSet::const_iterator c = components_.begin();
	for( size_t i=1; i <= size; ++i, ++c)
	{
		chemsys.X[i] = (*c).C;
	}
}

void Solution:: load( const array<Real> &C )
{
	if( C.size() < size )
	{
		throw exception("Solution::load(invalid array)");
	}
	ComponentSet::iterator c = components_.begin();
	for( size_t i=1; i <= size; ++i, ++c)
	{
		(*c).C  = C[i];
	}
	
}

void Solution:: save( array<Real> &C) const
{
	if( C.size() < size )
	{
		throw exception("Solution::save(invalid array)");
	}
	ComponentSet::const_iterator c = components_.begin();
	for( size_t i=1; i <= size; ++i, ++c)
	{
		C[i] = (*c).C;
	}
}

void Solution:: display( std::ostream &fp ) const
{
	fp << "/---------------" << std::endl;
	for( ComponentSet::const_iterator c = components_.begin(); c != components_.end(); ++c )
	{
		const Component &comp = *c;
		fp << "| [" << comp.species.name << "]=" << comp.C << std::endl;
	}
	fp << "\\---------------" << std::endl;
}

Real Solution:: sumCharge() const throw()
{
	Real sum = 0;
	for( ComponentSet::const_iterator c = components_.begin(); c != components_.end(); ++c )
	{
		const Component &comp = *c;
		sum += comp.C * comp.species.Z;
	}	
	return sum;
}


void Solution:: mul( Real a ) throw()
{
	for( ComponentSet::iterator c = components_.begin(); c != components_.end(); ++c )
	{
		(*c).C *= a;
	}
}


Real Solution:: pH() const
{
	return -Log10( (*this)[ "H+" ] );
}

void Solution:: zero() throw()
{
	for( ComponentSet::iterator c = components_.begin(); c != components_.end(); ++c )
	{
		(*c).C = 0;
	}
}

void Solution:: addLeak( const Real zeta_m, const Real CellVolume, const Solution &in, const Solution &out )
{
	assert(&in.library  == &library);
	assert(&out.library == &library);

	ComponentSet::iterator        iter_self = components_.begin();
	ComponentSet::const_iterator  iter_in   = in.components_.begin();
	ComponentSet::const_iterator  iter_out  = out.components_.begin();

	for( size_t i=size;i>0;--i, ++iter_self, ++iter_in, ++iter_out)
	{
		Component       &c_self  = *iter_self;
		const Component &c_in    = *iter_in;
		const Component &c_out   = *iter_out;
		const Species   &species =  c_self.species;
		
		assert( &c_self.species == &c_in.species);
		assert( &c_self.species == &c_out.species);
		
		const Real zzm  = zeta_m * species.Z;
		const Real Psim = Psi( zzm );
		const Real SPm  = species.SP(zeta_m);
		//std::cerr << "P_{" << species.name << "}(" << zeta_m << ")=" << SPm << ", Psi=" << Psim << std::endl;
		//std::cerr << "\tout=" << c_out.C << ", in=" << c_in.C << std::endl;
		c_self.C +=  species.Pf * Psim * SPm * ( c_out.C - c_in.C * Exp( zzm ) ) / CellVolume;
	}
	
}

Real Solution::  set_dLeak_dzeta( const Real zeta_m, const Real CellVolume, const Solution &in, const Solution &out )
{
	derivative<Real> drvs;
	ComponentSet::iterator        iter_self = components_.begin();
	ComponentSet::const_iterator  iter_in   = in.components_.begin();
	ComponentSet::const_iterator  iter_out  = out.components_.begin();

	Real ans = 0;
	for( size_t i=size;i>0;--i, ++iter_self, ++iter_in, ++iter_out)
	{
		Component       &c_self  = *iter_self;
		const Component &c_in    = *iter_in;
		const Component &c_out   = *iter_out;
		const Species   &species =  c_self.species;
		numeric<Real>::function SPP( &species, & Species:: get_SP_times_Psi );
		
		const Real zzm  =  zeta_m * species.Z;
		const Real Psim =  Psi( zzm );
		const Real SPm  =  species.SP(zeta_m);
		
		const Real part1 = drvs( SPP, zeta_m, 0.01 ) * ( c_out.C - c_in.C * Exp( zzm ) );
		const Real part2 = SPm * Psim * ( -species.Z * c_in.C * Exp( zzm ) );
		c_self.C = species.Pf * (part1+part2)/CellVolume;
		
		ans += species.Z * c_self.C;
	}
	
	return ans;
}



#include "yocto/math/kernel/algebra.hpp"

//! this is "slow" dS, compute new dS
void ChemSys:: updateRate( Solution &dS, Real t , const Solution &S)
{
	//--------------------------------------------------------------------------
	//-- use the current solution to compute matrices
	//--------------------------------------------------------------------------
	S.save_into( *this );
	//computeGammaAndPhi(t,true);
	computeGammaAndPhi(t,false);
	computeXi();
	computeChi();
	
	//--------------------------------------------------------------------------
	//-- load perturbation into X
	//--------------------------------------------------------------------------
	dS.save_into( *this );

	//--------------------------------------------------------------------------
	//-- time independent dampening
	//--------------------------------------------------------------------------
	algebra<Real>::mul( dX, Chi, X );   //-- dX = Chi * dS
	
	//--------------------------------------------------------------------------
	//-- time dependant dampening
	//--------------------------------------------------------------------------
	algebra<Real>::muladd(dX,Xi,dtGam); //-- dX += Xi * dGamma/dt

	//--------------------------------------------------------------------------
	// update perturbation
	//--------------------------------------------------------------------------
	for( size_t i=M; i>0; --i )
	{
		X[i] -= dX[i];                 //-- dS -= (Chi*dS+Xi*dGamma/dt)
	}
	
	//--------------------------------------------------------------------------
	//-- store modified X
	//--------------------------------------------------------------------------
	dS.load_from( *this );
}


void Solution:: write_header( ios::ostream &fp ) const
{
	for( ComponentSet::const_iterator c = components_.begin(); c != components_.end(); ++c )
	{
		const Component &comp = *c;
		fp(" '%s'", comp.species.name.c_str() );
	}	
}

void Solution:: write_values( ios::ostream &fp, const char *fmt) const
{
	assert(fmt!=NULL);
	for( ComponentSet::const_iterator c = components_.begin(); c != components_.end(); ++c )
	{
		const Component &comp = *c;
		fp(fmt, comp.C );
	}	
}

